Theory of STM junctions for 7r-conjugated molecules on thin insulating films 



Sandra Sobczyk, Andrea Donarini,* and Milcna Grifoni 

Institut fur Theoretische Physik, Universitat Regensburg, 93040 Regensburg, Germany 

(Dated: May 23, 2012) 

A microscopic theory of the transport in a scanning tunnelling microscope (STM) set-up is in- 
troduced for 7r-conjugated molecules on insulating films, based on the density matrix formalism. 
A key role is played in the theory by the energy dependent tunnelling rates which account for the 
coupling of the molecule to the tip and to the substrate. In particular, we analyze how the geo- 
metrical differences between the localized tip and extended substrate are encoded in the tunnelling 
rate and influence the transport characteristics. Finally, using benzene as an example of a planar, 
rotationally symmetric molecule, we calculate the STM current voltage characteristics and current 
maps and analyze them in terms of few relevant angular momentum channels. 
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I. INTRODUCTION 

Scanning tunnelling microscopy (STM) is an impor- 
tant tool for imaging surface structures and for studying 
the electronic properties of individual molecules since its 
introduction by Binnig and Rohrer ' . Various authors 
have developed theories of STM 3 ~ 8 ' 10 ~ 15 , among those 
the famous ones published by Tersoff and Hamann 4 ' 8 ' 9 
in the 1980s. Their work is the basic theory used to ex- 
plain STM images without atomic resolution 16 , i.e. STM 
images with characteristic feature sizes of > lnm, for ex- 
ample the scattered waves of surface states, as well as ad- 
sorbatcs, defects and substitution atoms on the surface 17 . 
Tersoff and Hamann showed that those experiments, as 
those on reconstructed Au surfaces, may have a simple 
explanation. In their articles the tip was modeled as a 
spherical potential well of radius R = 9 A, taking the s- 
wave solution of the macroscopic Schrddinger equation to 
describe the electronic tip-state. With Bardeen's pertur- 
bation theory of tunnelling 18 , they showed that the STM 
image is approximately the Fermi-level local density of 
states (LDOS) contour of the sample at the center of 
the sphere. Though the Tersoff-Hamann approach can- 
not be used to explain famous STM experiments that 
show atomic resolution, because it ignores the detailed 
structure of the tip wave functions. For true atomic res- 
olution, for which the length scale is much smaller than 
one nanometer, the convolution of tip states and sam- 
ple states must be taken into account 19 . Chen presented 
an extension of the Tersoff-Hamann theory that implies 
more detailed tip-models and allows to interpret higher 
resolution STM images 11 ' 20,21 . Several other authors sug- 
gested that atomic resolution demands small tip-sample 
distances 6 ' 10,22 , which are not fully described within the 
Bardcen tunnelling theory 18 . 

In fact the majority of the STM studies of single 
molecules, in experiment and in theory, has so far been 
limited to molecules on metals or semiconductors. In 
these cases the electronic properties of an individual 
molecule are strongly perturbed by the presence of the 
substrate electrons. In order to understand the electronic 
properties of an individual molecule, an electronic decou- 



pling from the supporting substrate is desirable. Hence, 
in the seminal experiments [23,24], STM measurements 
have been performed on molecules on insulating films 
having a thickness of only few atomic layers. The layer 
is in turn grown on top of a metallic substrate. This set- 
up allows to electronically decouple the molecule from 
the metallic surface, so that electronic properties of in- 
dividual molecules can be studied. At the same time the 
electrons can still tunnel through the insulating films, fa- 
cilitating imaging with the low-temperature STM at a 
low tunnelling current. 

In this work we present an STM theory that enables to 
study the transport properties of individual 7r-conjugated 
molecules in the latter STM configuration. We model 
the device with a double-barrier tunnelling set-up, and 
treat its dynamics in the sequential tunnelling limit via 
a density matrix approach. We show that the geometri- 
cal aspects in the coupling to the substrate and the tip, 
results into significantly different, energy dependent tun- 
nelling rates. Using benzene as an example, we calculate 
current voltage characteristics and constant height cur- 
rent maps for different biases and substrate work func- 
tions, thus simulating STM images with atomic resolu- 
tion. Due to the rotational symmetry of the benzene 
molecule we express the theory in the angular momen- 
tum basis, and we prove that the tunnelling dynamics 
from/to the extended substrate is described by angular 
momentum channels. Vice versa, the localized tip mixes, 
in the tunnelling events, the angular momentum states 
of the molecule. This mixing produces, for specific sub- 
strate work functions, negative differential conductance 
and current blocking also detectable in the topography 
of the STM surface plots. 

Both the Pauli and the generalized master equation 
have been repeatedly used in the modelling of STM 
junctions 25-30 . Nevertheless, to our knowledge, STM 
junctions with a thin insulating layer have not been sys- 
tematically studied within the framework of the general- 
ized master equation. 

This paper is outlined as follows: in section II we 
present a general transport theory for 7r-conjugated 
molecules in the STM set-up. We introduce the model 



Hamiltonian of the system and provide a detailed analysis 
of the tunnelling dynamics in terms of energy dependent 
tunnelling rates. In section III we apply the theory to 
a benzene molecule. The corresponding current voltage 
characteristics and current maps are discussed in section 
IV. Finally, conclusions and remarks are presented in sec- 
tion V. 
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II. LOW ENERGY THEORY OF STM ON 
INSULATING LAYERS 

A. Hamiltonian and tunnelling amplitudes 

A scanning tunnelling microscopy (STM) set-up with 
a thin insulating film involves the STM tip, the substrate 
and the molecule (Fig. la), weakly coupled to each other. 
Therefore we can describe the whole system by the total 
Hamiltonian 



H = H„ 



-Htip + H t 



(1) 



The first term gives the Hamiltonian of an arbitrary tt- 
conjugated molecule. We assume that only the 7r-orbitals 
contribute to transport. Thus, to each atom is assigned 
only one orbital (the 2p z orbital orthogonal to the plane 
of the molecule), while the entire a backbone is included 
only via the parametrization of the Hamiltonian for the 
7r-conjugated electrons. The latter, written in the atomic 
basis, is a simplified version of the Pariser-Parr-Poplc 
(PPP) Hamiltonian 31 ' 32 , expressed in terms of the non- 
interacting Huckcl-Hamiltonian 33 and a constant inter- 
action term: 

H m = agS^dga + b a pdt i aa dp<j + 

aa a=£fia 

+ \u{N-N f, 

where dJ a<7 creates an electron of spin a in the p z -orbital 
of the atom a, and a = 1, M runs over the M atoms 
of the molecule. The hopping energies b a p are assigned 
using the Slater-Kostcr method 34 with atomic param- 
eters and geometrical configurations obtained from the 
literature. The on-site energy for the atom a is de- 
noted by a a and can also vary from atom to atom. Fi- 
nally, the constant interaction model 35 assumes that the 
Coulomb interaction between the electrons is parameter- 
ized by a constant capacitance C, what is finally dcfin- 

2 

ing the Coulomb interaction U = where e is the 
charge quantum. This model also assumes that the dis- 
crete single-particle energy spectrum is unaffected by the 
interactions. Finally, N = d^ ac7 d aa - counts the num- 
ber of 7r-electrons in the molecule which is No for the 
neutral case. 

The simplicity of the Hamiltonian for the molecule pre- 
sented here allows to carry out most of the calculations 
(specifically the ones relative to benzene presented in sec- 
tions III and IV) at an analytic level since the many-body 
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FIG. 1: (Color online) Panel (a) - Sketch of the investigated 
STM set-up. A 7r-conjugated molecule, here benzene, is sep- 
arated by a metal substrate (yellow) only through an ultra- 
thin insulating film (red). A bias voltage is applied between 
the substrate and the tip. Panel (b) - Schematic illustration 
for the sum of the potentials of the substrate, the molecule 
and the tip v = v su t + v m + vu p , along the z direction. We 
choose the energy of the vacuum between the molecule and 
the tip, as well as the energy of the tunnelling barrier be- 
tween molecule and substrate to be zero. The energies at 
the bottom of the conduction band of tip and substrate are 
£ S/t _ _ < £S/t _ e y T ^ wnere £ S J T are the Fermi energies mea- 

S/T 

sured from the band bottom and $ are the work functions 
for the tip and the substrate. The work functions are shifted 
by the applied bias voltage. 



eigenstates of the interacting Hamiltonian coincide, in 
this case, with the ones of the non interacting one. Nev- 
ertheless, the transport theory is not affected by the par- 
ticular choice of the Hamiltonian for the molecule and the 
transport characteristics remain qualitatively the same 
for the different models, as far as the symmetry of the 
states is preserved. 

We consider the tip and and the substrate as reservoirs 
of non interacting electrons. In particular, we describe 
the metallic substrate as a potential well (see Fig. lb) 
with no confinement in the x and y direction. The asso- 
ciated Hamiltonian H su \, reads 



H. 



ub 



E4 



SJ 



C Ska C Skcr> 



(3) 



ha 



where = 



h 2 \k\ 2 



with k 



k„, k z ) and c'_- 



t 

Ska 



creates an electron of momentum k and spin a in the 
substrate and \zq\ is the z extension of the substrate (see 
Fig. 1). The continuous choice also for the z component 
of the momentum is justified in the limit |z | 3> where 
Xf is the Fermi wave length of the substrate. Only bound 
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states (e§ < 0) are considered in the calculation and their 

v k 

explicit wave function is given in the Appendix A. 

An analogous shallow square potential for the z direc- 
tion describes the metallic tip. A parabolic confinement 
in the x and y direction is though added to the model 
to simulate the spatial localization of the tip states. The 
tip Hamiltonian reads: 

ff tiP = Yl £ L C Tk z * C Tk;* , (4) 
k z a 

where = + Huj + and c^ Tk a creates an electron 
with momentum k z , spin a, and in the ground state with 
respect to the lateral confinement. 

We are confident that the particular choice of the con- 
finement for the tip Hamiltonian is not crucial for the 
results. Nevertheless, as it has already been theoretically 
predicted 11 and experimentally observed 36 , the symme- 
try of the tip is very important. We will restrict in this 
work to tip wave functions which are rotationally invari- 
ant with respect to an axis perpendicular to the surface 
of the substrate. 

The last term of Eq. (1) is the tunnelling Hamiltonian. 
It contains two parts: one for the substrate- molecule tun- 
nelling, the other for the tip-molecule tunnelling: 

\kia 

The index i denotes the molecular orbital, i.e. the lin- 
ear combination of the atomic p z orbitals introduced in 
Eq. (2), x = S, T indicates the substrate or the tip and we 
have introduced the general label k indicating the orbital 
quantum numbers of both the leads with the identifica- 
tion k = k for the substrate and k = k z for the tip. The 
coefficient €L is the tunnelling amplitude that contains 
all the geometrical information about the tunnelling pro- 

2 

cesses. Denoting by h = ^- + v m + v su \, + u t i P the single 
particle Hamiltonian for an electron in the STM set-up, 
we define this amplitude by 

i* := ( X ka\h\ia) , (6) 

where Ixfccr) and \ia) are eigenstates of the reservoir \ 
and of the molecule, respectively. The kinetic energy 

2 

of the electron is given by . The molecule, tip and 
substrate potentials are denoted by v m , Vti p and i> su b, 
respectively. The z-dependence of the total potential 
v = v m + w su b + Vtip is schematically shown in Fig. lb. 
It is the sum of three potential wells, for the substrate, 
molecule and tip where £q < defines the bottom of the 
conduction band and Eq + Ep < are the Fermi energies. 
For the tunnelling amplitudes, it follows: 

P 2 

f ki = (xk<j\ h v m + (xkv\v su b + v tip \ia) 

Zm s * 

~° (7) 
= Ei(xkcr\ia) = £i^^(xka\aa}(aa\ia) , 

a 



where h mo \ is the non-interacting single-particle Hiickel- 
Hamiltonian that satisfies the eigenvalue equation 
K w \\iv) = Ei\ia). 

The key observation to understand why the matrix ele- 
ment (xk<r\v SVL b + Vti p \ia) can be neglected while the con- 
tribution (xka\v m \ia) containing the molecular potential 
should be retained is the larger penetration length of the 
lead wave function, with respect to that of the molec- 
ular orbital, into the barrier region separating the lead 
and the molecule. This difference implies in fact that 
the relevant integration region for the matrix element 
(x/cer|v su b + v t ip + w m|*cr) is shifted towards the molecule. 
Consequently the kinetic energy contribution should be 
complemented by the one of the molecular potential. For 
systems characterized by states with comparable pene- 
tration lengths instead, the relevant integration region is 
in the tunnelling barrier and the kinetic energy yields the 
dominant contribution. 

The different penetration lengths for the lead and 
molecule wave functions is justified as follows. First, the 
spatial extension of the valence orbitals is larger for the 
metallic atoms of the lead than for the ones in the con- 
jugated molecule. Moreover, the states in the lead which 
dominate the tunnelling have no nodal planes perpendic- 
ular to the molecular surface (low fcy) while the HOMO 
and LUMO states of a conjugated molecule have usu- 
ally several nodal planes perpendicular to the plane of 
the molecule. These perpendicular nodal planes are as- 
sociated to a destructive interference between the atomic 
wave functions which implies that the higher the num- 
ber of nodal planes, the shorter is the extension of the 
molecular orbital in the direction perpendicular to the 
molecular plane. 

Notice that the energy of the vacuum between the 
molecule and the tip has been set to zero. Likewise we 
also set to zero the top of the tunnelling barrier between 
the molecule and substrate, corresponding to the thin 
insulating layer. The theory is not affected though by a 
different value of the potential in the barrier regions as 
far as the latter is spatially uniform if compared with the 
product of the lead and molecule wave functions in the 
same region. A more precise description of the lead po- 
tential would in first approximation just lead to a renor- 
malization in Eq. (7) of the orbital energy e,. 

In the last step of Eq. (7) we added the completeness 
1 = ^ a |aa)(acr|, where \ota) is the p 2 -state of the atom 
a, thus showing that the wanted matrix element can be 
expressed in terms of the overlap (xka\aa) of the lead 
and the p 2 -orbital and the basis transformation (aa\ia) 
from the molecular to the atomic orbital. Finally, we 
obtain for the tunnelling amplitudes: 

tki = e*X! 0r ^ 2 '^ ti P - Ra){ a(J \ ia ) , ( 8 ) 



t S ki = e i Y J ^ Aa O s {k){aa\icx) , (9) 

a 
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TABLE I: Parameters at, /3i used for the Gaussian p z -orbitals 



i 


i 


2 


3 




0.368 


1.113 


4.997 




0.502 


1.438 


2.620 



where R a and i?ti P are the position of the atom a and 
of the tip, respectively. The overlaps O x are given ex- 
plicitly in the Appendix A and are calculated using the 
p z -orbital 37 : 

p z (r-R a ) = (?\a G ) =n G Y / fo(r-R a )-e z e~ a '\ p -^ 2 , 

(10) 

where e z is the versor in the direction perpendicular to 
the molecular plane, the coefficient n G assures normaliza- 
tion and the parameters on and /3j , that we show in table I 
for the specific case of a carbon atom, define the gaussian 
representation for a Slater type orbital commonly used in 
DFT calculations 38 - 39 . Analo gous parametrizations are 
available also for other atoms and allow a straightforward 
application of the model to generic planar 7r-conjugated 
molecules. The overlap functions of the substrate and 
the tip are qualitatively different since they reflect the 
different geometries of the corresponding contacts. The 
plane wave description of the electrons in the substrate 
implies that in Eq. (9) the position of the atom R a only 
appears in the phase factor as a scalar product with the 
component of the momentum parallel to the substrate, 
fc| | . Additionally we obtain a function that only depends 

on the electron's momentum k in the substrate and on the 
thickness of the insulating barrier. This particular form 
already suggests that the tunnelling between the sub- 
strate and the molecule is not an incoherent collection of 
tunnelling events happening in correspondence to the dif- 
ferent atoms since their position is recorded in the phase 
of the tunnelling amplitude. Some of the consequences 
of this spatial coherence will appear more clearly in sec- 
tion III where we analyze the special case of a benzene 
STM junction. The overlap function for the tip is more 
complex. Due to the cylindrical symmetry of the tip and 
atomic orbital with respect of their rotational axes, we 
can only further conclude that only the modulus of the 
component of i?ti P — R a parallel to the molecular plane 
influences the tunnelling (see Appendix A). 



B. Tunneling dynamics 

Our method of choice to treat the dynamics in the 
regime of weak coupling between system and leads is the 
Liouville equation method. A detailed discussion and 
derivation of the equations of motion for the reduced den- 
sity operator of the system can be found e.g. in 40,41 ; we 
will give here only a short overview adapted to the STM 
set-up. 

We start from the Liouville equation for the total den- 
sity operator p(t) of the whole system consisting of the 
molecule, the tip and the substrate. Using the interac- 
tion picture and treating the tunnelling Hamiltonian (5) 
as a perturbation we get: 

in ^T = Kn(t),P I (t)}, (11) 

where the subscript I indicates the use of the interaction 
picture. Since we are not interested in the microscopic 
state of the leads, we focus on the time evolution of the 
reduced density matrix (RDM) a = Trs+r{p(t)}, which 
is formally obtained by taking the trace over the unob- 
served degrees of freedom of the tip and the substrate. 
The equation of motion for the RDM reads to lowest non- 
vanishing order in the coupling to the substrate and the 
tip 42 

a = ~[H m ,a] - -[H eS ,a] + C tun a := £<r. (12) 
n h 

The first term of this so called generalized master equa- 
tion (GME) gives the coherent evolution of the system in 
absence of the substrate and the tip. In the secular ap- 
proximation we only keep coherences between degenerate 
states and thus this term vanishes 40 . The commutator 
with H c f[ includes the normalization of the coherent dy- 
namics introduced by the couplings to the leads. Finally, 
the operator £tun describes the sequential tunnelling pro- 
cesses. The sum of these three contributions defines the 
Liouville operator C. 

Let us concentrate first on the tunnelling processes oc- 
curring in the system. The corresponding contribution 
to the master equation, projected into the subspace of 
A^-particles and energy E reads: 



NE 



V 



NE 



4r r U E - H ™)fx(E - H m )d JT + d JT T?AH m - E)f+(H m - E)d 



NE 



E E [4^ - E')a N -™'f+(E E')d jT + d jT Y*{E' E)a N+1E ' f~ (E' E)d\ T 

XT ijE' 



Vne 
(13) 



where a := VnevVne being Vne '■= J2i \NEl)(NEl\ the projection operator on the subspace of N particles 
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and energy E, and I the additional quantum number 
that distinguishes between degenerate states. Moreover, 
is the Fermi function for the lead x, fx( x ) := 
f(x — fi x ), and fv( x ) '■= 1 — fx( x )- The terms pro- 
portional to f£{x) describe in (13) tunnelling events to 
the molecule, while the tunnelling out of the molecule is 
associated to fz{x). Finally fi x stands for the electro- 
chemical potentials of the substrate or the tip. They are 
defined via the applied bias voltage as fi$ = Mo + (1 — 
c)eVb, ht — Mo — c eVb and consequently eVb = /is — Ht , 
with the electron charge e, the equilibrium potential [j,q 
and the coefficient c governing the relative bias drop 
at the tip and the substrate. A symmetrical potential 
drop is obtained for c = 1/2, while for c = 1 the bias 
drops completely at the tip-molecule interface. Finally 
/in = —&o relates the equilibrium chemical potential to 
the work function and, in equilibrium, the work func- 
tions of the two leads are assumed equal. Beside the 
Fermi function, the tunnelling rates are characterized by 
the geometrical component: 



(14) 



The argument AE of the rate Tfj is the energy difference 
En+i — En of the many body states involved in the tun- 
nelling process, sometimes written in Eq. (13) in terms 
of the operator H m . Notice that the rate vanishes if 
AE > since we restrict the Hilbert space of the leads 
to the bound states i.e. Sk < 0. The quantity r* plays a 
central role in the theory and in the following section we 
will discuss its calculation in detail for the tip and the 
substrate case using the example of a benzene molecule. 

A natural expression for the current operators is ob- 
tained in terms of the time derivative of the reduced den- 
sity matrix: 



(/sub + /tip) =J2 Tr i N * NE ] . 



(15) 



NE 



where L, 



.b/tip 



are the current operators calculated for 



the substrate and the tip interfaces. Conventionally we 
assume the current to be positive when it increases the 
charge on the molecule. Thus, in the stationary limit, 
(/sub + /tip) is zero. The stationary current is obtained 
as the average: 

(7 sub ) = Tr{cr sta t/ S ub} = -(/tip) , (16) 

where er s t a t = limt-s-oo a(t) is the stationary density op- 
erator that can be found from 



Cstat 



0, 



(17) 



where C is the Liouville operator. Finally, by following 
exactly the procedure given in [41], we find the explicit 
expressions for the current operators: 



h= E Vne 

NEaij 



d j(T Tf j (H in -E)f+(H m -E)dt 



dlTf 3 (E-H m )f;(E-H m )d 



J a 



(18) 



Vne 



where the energy renormalization terms, present in the 
GME, do not appear. 

Since the tunnelling changes the number of electrons 
on the molecule, the latter behaves as an open system and 
it is useful to introduce the operator H' m = H m — /j,qN 
where N counts the number of electrons on the molecule. 
For example, at zero temperature and zero bias the equi- 
librium is reached when the molecule is in the ground 
state of H' m and not of H m . As we have already 
shown elsewhere 43 , also the non-equilibrium conditions 
for transport can be better understood in terms of the 
spectrum of H' m . For this reason in Figs. 3 and 4 the 
geometrical part of the rates is plotted as a function of 
AE' := AE- mo- 



III. THEORY APPLIED TO BENZENE 

The molecular orbitals of benzene are also eigenfunc- 
tions of the projection / of the angular momentum along 
the main rotational axis, which we assume to be the z- 
axis. Therefore, the basis transformation that occurs in 
Eq. (7) reads for a benzene molecule 



(aa\la) 



1 

7T 



(19) 



and the corresponding single particle eigenenergies E\ , oc- 
curring in the Eqs. (8) and (9) for the tunnelling ampli- 
tudes, read: 



, 2?r, 

Ei = a + 2b cos ( — I 



(20) 



For a benzene molecule the possible values of the an- 
gular momentum quantum number I are 0,±1,±2,3 
corresponding to the energy level scheme of the Hiickcl 
Hamiltonian shown in Fig. 2. Since the Hamiltonian 
is invariant under the discrete rotations of angles ?i7r/3 
with n € Z, the same quantum numbers also label the 
many-body eigenstates of the benzene molecule, irrespec- 
tive of the complexity of the description of the Coulomb 
interaction 41 . All the single particle states show a twofold 
spin degeneracy but only few states possess an additional 
twofold orbital degeneracy. The latter is essential for the 
explanation of the transport features of benzene within 
an STM experiment. 



A. The substrate-molecule tunnelling rates 

Let us start with a detailed discussion of the substrate- 
molecule tunnelling rate. To perform the sum over the 
momenta k in Eq. (14) we transform it into energy inte- 
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FIG. 2: Energy levels of the Hiickel Hamiltonian and the 
corresponding values of the angular momentum 



grals, using the definitions en := h ^ ^ and e z := 



E = EE> 



(21) 
(22) 
(23) 



where the volume V = \zq\S is canceled out in the ther- 
modynamic limit by the normalization of the orbitals 
which define the overlap function. Moreover we observe 
that Eq. (14) requires the calculation of the product 



fczy w 



si S i>J2\ s(k)\ 2 (aa\la}(a'a\l'a) 



,+ifc, i 



(24) 



an' 



We write the exponential function in Eq. (24) as 
e +i%„-(R a -R a ,) = e +i\k u \\R a -R a ,\costi and thc cquation n _ 

nally becomes 



(25) 



x e 



+i V /^£|||Afl 7 |co S1 ? 



|Os(e N ,ez)| : 



where we introduced a — a' := 7 , |i2 Q — = |Ai? 7 |. 
We insert Eq. (25) in the substrate case of Eq. (14) and, 
after solving the integral over d?9, we find: 



x(l-l') 



f £F+*0 /•£? + *<] y / j2 m 

x / d£|| y„ ^^E^iv^iii^i 



x |Os(ei|,e,)| 2 e+ i * , 'T% Jt -A£;) 



(26) 



with Jo(x) the zero-order Bessel function. Finally, using 
the relation 



£ e± i^a(i-0 =6( j M/) 



(27) 



and the fact that e 1 ^ ' 7 = ^ 7 e 1 a «T ll the integral 
over Em yields 



*E^(7¥(^ 

7 \ 



r 



e,-4')|Ai? 7 | Je-'T'T (28) 



x |O s (A£;- e2 -4, ez )| 2 
x6(AB-e z - eg) 6 (e 2 - A_E) . 

The integral in Eq. (28) has to be solved numerically. 
The main result of the latter calculations is 



SwTf(AE) 



(29) 



which ensures that tunnelling processes involving the 
substrate are happening through angular momentum 
channels because a mixing of angular momenta is not al- 
lowed in the substrate. We will see that this only happens 
for substrate-tunnelling-processes, while there is no con- 
servation rule for angular momenta in the tip-tunnelling 
case. The function Tf(AE) is the geometrical rate and 
we plot it in Fig. 3 for different angular momenta. The 
rates decrease of several order of magnitudes by increas- 
ing the absolute value of the projection of the angular 
momentum I. This is the direct consequence of the de- 
creasing extension of the molecular orbitals in the direc- 
tion perpendicular to the molecular plane with increasing 
the number of vertical nodal planes. 

The lower limit of the energy axis in Fig. 3 is 
the upper limit is the work function c/>q . These limits are 
set by the substrate model in which only bound states 
of a single band are taken into account (ejf < ef < 0). 
While approaching the low energy limit AE = —e§ both 
the density and the penetration length of the states in 
the substrate which contribute to the rate reduce, hence 
the turn down. On the other hand, the increasing of the 
density of states and of the penetration length explains 
the turn up at the upper energy border (AE = (p$ ). 



e§ while 



B. The tip- molecule tunnelling rates 

Let us now discuss the tunnelling events happening be- 
tween the tip and the molecule. To model the tip we con- 
sider a harmonic confinement in the x and y directions. 
By considering the tip to be in the ground state of the 2- 
dimensional harmonic oscillator, the longitudinal energy 
Em is fixed to be thc constant em = tuv, cf. below Eq. (4). 
The sum in Eq. (14) thus transforms into a sum over 

k z . Because of the relation k z = \I^S-e z we can replace 
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-6 -4 -2 2 4 
AE' = AE-|J [eV] 

FIG. 3: (Color online) Tunneling rate Tf describing 
substrate-molecule tunnelling processes for different angular 
momentum quantum numbers I. The thickness of the sub- 
strate barrier is d = 3A, while work function and Fermi energy 
are respectively (j>o = AeV and ef = 7eV 



the sum by the integral: £) fc -> j-^/f J de z |Zcnd Ztipl 



Eq. (8) implies 



x 0^(k z ,R tip - R a )0 T {k z ,Rti P - R a > 



(30) 



that we insert in Eq. (14). After solving the energy inte- 
gral we finally find 



Tf l ,(AE,R tip ) = w J-Y^e mi i 



x 0* T (k, Ra p - R a )0 T (k, R tip - R a .) 



i^g-{al~a'l') 

I ^end ^tip | 



^/AE -tujj 



x Q(AE -fiu- el)Q{2hw - AE + ej) , 

(31) 

where k = j^(AE — tuo — e^). The occurrence of 

both I and V in the latter equation, shows that a mix- 
ing of angular momenta during the tip-tunnelling process 
takes place. Upon inspection of Eq. (31) we find some 
important relations obeyed by the tunnelling rate, where 
we use the fact that I and V always occur in the form 
V = ±1: 



rl = ?l = (if)* e 



(32) 



where Tf = T^, which implies the existence of an angu- 
lar momentum dependent phase when I =/= V . In Fig. 4 
we show the diagonal elements of the rate matrix Tj v ex- 
emplified for I = ±1 and I = ±2. As for the substrate, 
the channel I = ±1 leads to a much larger rate than the 
channel I = ±2. The phase in the off diagonal elements 
depends on the tip position i? t i P and it is calculated as 



(-^tip) 



arg(ifc-)- 



(34) 



In Fig. 5 we show the values acquired by the phase 
<M-Rtip) a s a function of the tip position. The phase is 
approximately constant along the radii leaving the cen- 
ter of the molecule. Due to the cylindrical symmetry of 
the tip wave function a good approximation to the phase 
(f>l(Rtip) is given by: 



<M-Rtip) = Wtip, 



(35) 



where #tip is the angle describing the projection of the 
tip position on the molecular plane if the origin is the 
center of the molecule. By convention we assume # t ip = 
along the radius that intersects the position of the atom 
of the molecule (see Fig. 5). The derivation of this 
simple expression for (pi as well as a discussion on its 
limits of validity are given in Appendix C. Notice that 
the phase defined in Eq. (34) only depends on i? t ip even 
if O a contains the bias. Nevertheless, the tunnelling rate 
Eq. (33) depends on the bias via the Fermi energy. 

In Fig. 5 the position of the <pi = line is arbitrary 
and connected to the arbitrary choice of overall phase 
for the molecular orbital with angular momentum I. A 
different choice of the overall phase would, nevertheless 
simply appear as a rigid rotation of the plots. Moreover, 
this arbitrariness has no influence on the current voltage 
characteristics of the junction. 

In the substrate the tunnelling matrix is diagonal and 
proportional to the identity matrix, independent of the 
basis representation, see Eq. (29). In contrast, according 
to Eq. (33), off-diagonal elements are present in the tip- 
tunnelling matrix which, in the basis {|Z), \l)}, reads 



r T = rr 



1 e" 

e +2i<M_R tip ) 




(36) 



An interesting effect of the localized character of the 
tunnelling from/to the tip can be better appreciated by 
switching to the basis which diagonalizes the matrix in 
Eq. (36). The substrate rate matrix is still proportional 
to the identity matrix. For the tip rate matrix we get 
instead: 



where we have introduced the notation I = —I. Thanks 
to the relations (32) we can rewrite the tunnelling rate 
as 



(33) 



r T = rr 



2 0' 
.0 0, 



(37) 



One diagonal element becomes zero, indicating that 
there are states which are coupled to the substrate but 
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N 

X 




-2 -1 

AE' = AE - no [eV] 

FIG. 4: (Color online) Diagonal elements Ff of the tip tun- 
nelling rate matrix Yj v for the different angular momentum 
states. The rates are calculated assuming z t i P - d = 3.5 A, 
4>q = AeV , ep = 7eV and ftu = AeV. The presence of the 
harmonic confinement explains also the different energy lim- 
its with respect to the ones of Fig. 3. The lower limit is at 



late the time evolution of the reduced density matrix as- 
sociated to £tun and the corresponding stationary state. 
The stationary density matrix is block diagonal in parti- 
cle number, energy and spin. In particular, if we restrict 
the dynamics to low biases, the only relevant states en- 
tering the dynamics are the states \bglr), |6 g 00), and 
\7glr), being the cation, neutral and anion ground states 
respectively. The neutral ground state is non degener- 
ate while the anion and cation are four times degenerate, 
due to the combination of the spin and orbital degenera- 
cies. The specific form of the stationary density matrix 
depends on the bias, the temperature, and the tip posi- 
tion. Nevertheless, due to the form of the tunnelling rate 
matrices, the two dimensional sub-blocks corresponding 
to orbitally degenerate states have always the following 
structure: 



NE a 



A £ e -2i<MB tip ) 

£ e +2i<^(-Rtip) A 



(38) 



-ep + hu while the upper limit is at 



-ep + 2hu). 




x[A] 



x[A] 



FIG. 5: (Color online) Phase <f>i of the tunnelling rate matrix 
Yj v , Eq. (34). The phase is almost constant if the tip is moved 
along the radii outgoing from the center of the molecule. The 
carbon atoms are labelled by a = 0, . . . , 5. 



not to the tip. The decoupled state represents a blocking 
state, which can be populated by a tunnelling event from 
(to) the substrate but cannot be depopulated by a tun- 
nelling event to (from) the tip. The presence of blocking 
states is visible in the current-voltage characteristic, as 
we will discuss in the next section. 



where N = 5,7, the spin r =t, I and the parameters A, B 
are functions of the tip position i?ti P and of the bias Vb 
(see Appendix B). This result is a posteriori not surpris- 
ing. The comparison of Eq. (38) with Eq. (36) reveals 
in fact that the density matrix and the rate matrices are 
diagonalized by the same basis transformation (the sub- 
strate rate matrix is diagonal in all bases). Thus, the 
form of (7 s tat could be calculated from the observation 
that the dynamics of the populations and the coherences 
is decoupled when expressed in the eigenbasis of the rate 
matrices. It should be noticed that the diagonalizing ba- 
sis depends on the phase, see Eq. (34), which in turn 
depends on the tip position. Thus it is not possible to 
describe the system using only populations in a unique 
basis valid for all the positions of the tip. 



D. The effective Hamiltonian 



C. Stationary density matrix 

By combining now the expression for the tunnelling 
rates with the dynamical equation Eq. (13) we can calcu- 



Until now we only concentrated on the sequential tun- 
nelling processes in the system. We still have to discuss 
the imaginary term in Eq. (12) which contains the effec- 
tive Hamiltonian i? c fr- The latter is defined as: 



H ^ = 2^ J2 EE Vne \ d t T u> ( E - H m)Px ( E - H m )d Va + d Va Tl, (H m - E)p x (H m - E)d} la \ V NE , (39) 



NE x<? U' 



with the projector Vne = J2 n \NEri)(NEn\ 
and the principal part functions Px( x ) = 



Re* 



2 + 2ixk B T ( X Px) 



with T being the tern- 
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perature and VP the digamma function. Eq. (39) shows 
that the effective Hamiltonian is block diagonal in 
particle number and energy, exactly as the density 
matrix in the secular approximation. Consequently, it 
only influences the dynamics of the system in presence 
of degenerate states with corresponding subblocks larger 
than a mere complex number. For the sake of simplicity 
we will include in the following calculations only the an- 
ion ground states, (i.e. the spin and orbitally degenerate 
7 particle ground states). Analogous arguments holds 
for all the other degenerate states of the molecule. 

If Fiii cx 5w (substrate case, see Eq. (29)), the effective 
Hamiltonian iJ c ff in the 7 particle ground state subspace 
is proportional to the identity matrix, as can be proven 
from Eq. (39) remembering that H m conserves the an- 
gular momentum and it is invariant under the symmetry 
operation that brings \7 g lr) into \7 g lr) and moreover that 
= r|f. Thus, the substrate contribution to i? c ff triv- 
ially commutes with er s t a t ■ If the angular momenta I and 
I' can mix, like in the tip case, H e s acquires off diagonal 
terms and a more detailed discussion is required. In par- 
ticular, the form of the off diagonal elements depend on 
the particular model taken to describe the interaction on 
the molecule. As shown in the Appendix D, within the 
constant interaction model, the effective Hamiltonian for 
the tip can be written in the form: 



off 



(40) 



AE'5-6 




+2 1-2 



Particle Number 



where 



FIG. 6: (Color online) Together with a change in the en- 
ergy, the transition from the 6-particle ground state to the 
7-particle (5-particle) ground states is also associated with a 
change in the angular momentum of Al — ±2 (Al = ±1). 



IV. I-V CHARACTERISTICS AND CURRENT 
MAPS OF A BENZENE MOLECULE 



In the following discussion of the current voltage char- 
acteristics and current maps we only consider the ground 
state transition |6 S 00) O \7 g lr) or |6 9 00) «-> \h g lr). In 
Fig. 6 we represent the corresponding energy levels as a 
function of the particle number for a particular choice of 
the work function (we assume $ J = <I>q so that the chem- 
ical potentials are the same at Vb = 0). In the tunnelling 
event the molecule changes its particle number, angular 
momentum and energy (see Fig. 6). All these changes 
leave their fingerprints in the current voltage character- 
istics and current maps presented in Figs. 7-9. 



w = -(7>|4|6 s 0)(6 g 0|d r(T |7» 

7T 



xTf(E 7g -E 6g )p T {E 7g -E ( 



(></ ) 



+ -(7 g la\d llT \8 g 02a}(8 g 02a\dl\7 g la) 



Tf(E 8g -E 7g )p T (E 8g -E 



7g) 



(41) 



is the renormalization of the Bohr frequencies for the 
system and 



L = ±( 1 

2 I e +2i0 i (iJ ti p) 




(42) 



* n = 5.0 eV, T = 8K: 





Kja) V b = -1.688 V 




\ 

b) Coulomb- 


1 o 




blockade 




\ 


I w 




a) resonance 









V b [V] 



[pA] 



b)V b = -1.680 V 



o 



x[A] 



l[pA] 



FIG. 7: (Color online) Current voltage characteristics and 
current maps associated to the neutral-anion transition. The 
current maps are calculated with zti P — d = 5 A. Notice that 
the map in the Coulomb blockade region is just a rescaling of 
the one at resonance. 



Hence the effective Hamiltonian H^ s commutes with 
the stationary density operator cr s t a t given in Eq. (38). 
In conclusion, even if different from zero, the effective 
Hamiltonian does not contribute to the stationary dy- 
namics of our system because it commutes with the sta- 
tionary density matrix Eq. (38) calculated using only the 
tunnelling component of the Liouvillean. For a generic 
description of the Coulomb interaction on the molecule, 
corrections to H c g given by the 8 and 6 particle excited 
states should be taken into account and the form of H g 
is modified. For the sake of simplicity we restrict here 
to the constant interaction model. More details on the 
derivation and the discussion on the most general case 
are given instead in the Appendix D. 



In particular, the current is exponentially suppressed 
at small biases (the so called "in gap region" of trans- 
port) due to the Coulomb blockade 44 . The bias at which 
current starts to flow corresponds to a resonant condition 
between the chemical potential in the source (or drain) 
lead and the difference in the energy AE between the 
many-body states participating to the transport. For 
this reason the current voltage characteristics (and the 
associated differential conductance traces) recorded with 
an STM junction represent a valuable spectroscopic tool 
to investigate the many-body spectrum of the molecule. 
One has to keep in mind nevertheless that i) the reso- 
nant bias depends on the value of the work function of 
the leads, ii) the bias drops very asymmetrically at the tip 
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and substrate interface with an associated very different 
amount of energy available to the molecular transition. 
The shift in the position of the resonance with the work 
function can be observed by comparing the positions of 
the step in the current at negative biases in Fig. 7 and 9. 



* n = 7.0 eV, T = 8 K: 







b) resonance 

a) Coulomb- 
blockade 





a) V b = 2.150 V 



l[pA] 

j 0.025 

0.02 
■0.015 



l[pA] 



b) V b = 2.156 V 



VbM 



<[A] 



x[A] 



FIG. 8: (Color online) Current voltage characteristics and 
current maps associated to the neutral-cation transition. The 
current maps are calculated with zti P — d = 5 A. Notice that 
the value of the current at resonance is much higher than the 
one relative to the neutral-anion case (see Fig. 7). 

In Fig. 9 one can also observe how the same molecu- 
lar transition (between the neutral and anionic molecule) 
gives signals at different biases if triggered by a substrate 
(Vb > 0) or a tip (VJ, < 0) tunnelling event. A larger bias 
(in absolute value) is needed for a substrate transition 
since most of the bias drop concentrates at the tip in- 
terface. Moreover the current signal obtained at positive 
bias is a peak instead of a step due to an interference 
blocking effect analogous to the one discussed in [41]. In 
the interference blocking region the system is blocked into 
a particular linear combination of the 7 particles ground 
states that can be populated from the substrate but can- 
not be depopulated towards the tip. 

The angular momentum channel involved in the trans- 
port depends on the difference in the angular momentum 
of the many-body states participating to the tunnelling 
events. The neutral-anion and neutral-cation transitions 
correspond to Al = ±2 and Al = ±1 respectively, cf. 
Fig. 6, thus involving different angular momentum chan- 
nels. Since the lower is the angular momentum of the 
channel the larger are the rates, the current associated 
to the neutral-cation transition is larger than the one of 
the neutral-anion one, as it can be seen by comparing 
the resonant currents of Figs. 7 and 8. By comparing 
the same figures one finds also qualitative differences in 
the constant heights current maps: yet another finger- 
print of the different states involved in the transitions. 
The same differences are also confirmed by the constant 
current images presented in Fig. 10. 

Finally, the current maps presented in Fig. 9 suggests 
that also the interference effects have a topographic sig- 
nature. The current map taken in the Coulomb blockade 
region is in fact qualitatively different from the one taken 
in the interference blockade. 

To conclude, a comparison with the widely applied Ter- 
soff and Hamann (TH) theory 4 ' 8 ' 9 is compulsory. In par- 
ticular, for what concerns the current maps presented 
in Fig. 7 and Fig. 8, we do not expect qualitative dif- 



ferences between the effectively single particle TH the- 
ory and our many-body approach. Yet, this is almost 
accidental for the following reasons: i) we decided for 
simplicity to describe the system using a constant inter- 
action model in which the many-body states are single 
Slater determinants; ii) the initial and final many-body 
states of the tunneling event (e.g the neutral and anion 
ground states) fix the corresponding variation of angu- 
lar momentum (Al — ±2). Consequently, in the par- 
ticular case of benzene, only one single particle orbital 
contributes to the current formula given in Eq. (18). In 
general, though, many Slater determinants are necessary 
to identify a single many-body state and many molecu- 
lar orbitals would contribute to the transport. Moreover 
TH would not be able to address the interference block- 
ing regime and the associated current maps since it is 
effectively a non interacting single particle theory. 



ct> n = 4.0 eV, T = 8K: 



a) Coulomb- 
blockade 



b) interference- 
blockade 




3 3.05 3.1 3.15 3.2 



o 



l[pA] 
0.03 



x[A] 

I [pA] 

□ I 0.02 
0.015 
|0.005 



x[A] 



FIG. 9: (Color online) Current voltage characteristics and 
current maps associated to the neutral-anion transition. In- 
terestingly the current map in the interference blockade re- 
gion shows novel topographic features if compared with other 
maps involving the same states (see also Fig. 7) . In the inset 
a zoom on the interference current peak is presented. 





FIG. 10: (Color online) Constant current topographic images. 
The left panel refers to the neutral-cation resonance, (4>o = 
7eV, V b = 2.156V, I = 300pA ), the right panel, instead, 
to the neutral-anion resonance ((/>o = 5eV, Vb = —1.688V, 
/ = WOpA). 
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V. CONCLUSIONS 

In this paper we presented an STM transport theory 
sufficiently general to be applied to any device consisting 
of a 7r-conjugated molecule weakly coupled both to the 
substrate and the tip. While the weak tunnelling cou- 
pling to the tip is a natural assumption in STM exper- 
iments, the weak coupling to the substrate is motivated 
by recent STM set-ups with substrates covered by a thin 
insulating film 23,24 ' 36 . 

The essentially different geometry of the STM tip and 
the substrate is reflected in the respective tunnelling am- 
plitudes, whose energy dependence induces, within a den- 
sity matrix approach, characteristic non-constant tun- 
nelling rate matrices. The latter play a central role in 
the Liouville operator, which determines the dynamics 
of the system, and in the current operator. 

Interestingly, for these system, due to the different pen- 
etration lengths of the metallic states of the tip/substrate 
and the molecular orbitals into the corresponding tun- 
nelling barriers, the tunnelling amplitudes cannot be cal- 
culated using the standard Tcrsoff and Hamann approach 
and an alternative method is proposed. 

As an application of our general results we used a ben- 
zene molecule that enabled us to express the theory in 
the basis of the angular momentum I. The explicit cal- 
culation of the tunnelling rate matrices in the momen- 
tum basis shows a fundamental difference between the 
tip and substrate tunnelling dynamics. The delocalizcd 
tunnelling at the substrate happens via angular momen- 
tum channels (diagonal tunnelling matrices) while the 
localized tip mixes the angular momenta (off diagonal 
matrices). 

A direct consequence of this different tunnelling sce- 
nario for the two leads is found in the current voltage 
characteristics. At voltages sufficiently large to lift the 
Coulomb blockade, interference blocking occurs when de- 
generate states participate to the transport. While the 
presence of degenerate states is a necessary condition for 
the interference, only the tip tunnelling can detect it due 
to its localized nature which mixes the angular momenta 
in the tunnelling event. 

Moreover, also the STM surface-images can be calcu- 
lated within our theory. By varying the work function 
of the substrate we show simulations of STM constant 
height current maps and constant current topographic 
images in which the transport is dominated either by 



neutral-anion or neutral-cation transitions. In particular, 
striking is the difference in the current maps obtained in 
the resonant and interference blocking regime although 
the same many body states participate to the transport 
(see Figs. 7 and 9). 
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Appendix A: Calculation of the overlap functions 

To calculate the tunnelling amplitude in equation (7) 
we need to calculate the overlap between the metal's wave 
function and the p z -orbital. The latter is given, in the 
Gaussian description by 37-39 

(f\a G ) = n G J2 Pj (r~ R a ) ■ e z e~ a ^-R a \ 2 ; (A1) 

3 

where ng is the normalization factor which ensures 
J d?^ (r \a) | 2 = 1, R a is the position of the atom a and e z 
is the versor in the direction perpendicular to the plane of 
the molecule. Since the overlap is calculated as a function 
of the quantum number k defining the lead wave function, 
we will call the bracket (xka\aa) overlap function. 



FIG. 11: Scheme of a 1-dimensional, finite potential well with 
borders a and b and depth Vq. 

In our model both the tip and the substrate are de- 
scribed in the z direction as potential wells 45 . For future 
reference we report here the general expression for the 
cigenfunction of an arbitrary 1-dimensional potential well 
of depth Vq and whose borders are a and 6, see Fig. 11: 



c -ko [jj s [ n (k z a) + cos(A: z a)] e KZ , if — oo < z < a 
vl/ fc= (z; a, 6, Vq) = n z sm(k z z) + cos(k z z) , if a < z < b (A2) 

e +Kb ^jj sm (& z &) _|_ CO s(k z b)] e~ KZ , if b < z < oo , 



I 

where n z ensures the normalization J dz\^fk s ( z )\ 2 = 1 an d 

U 



k z sin(fc 2 &) — ncos(k z b) 
k z cos(fc z 6) + n sin(k z b) 



12 



The occurring wave number reads k z 



2m, 



and substrate overlap functions. The integral is: 



■j-rVo — kz, respectively. Due to the large size 

of the potential well compared to the Fermi wavelength 
we neglect the quantization of k z obtained by the corre- 
sponding eigenvalue equation. 

We conclude this introductory part with the explicit 
calculation of an integral common to both the tip and 



F kz (a, b, V ,atj) = / z^! k __ (z + d; a, b, ^e^ 2 , 

./ — OO 

(A3) 

where for simplicity we have omitted in F the dependence 
on the parameter d. The integration yields: 



F km {a,b,V ,aj) = — ^ 



4a \ 



x ^ e 2Rc 



-\h z d 



(1 + kV) 



-ctj ( a— d+7 



-a j lb-d+-£±- 



















— erf 





















+ e 4Q , 



-Ae 



erf 



^-d+^y _ K ^ 



Bc- Kd | 2^/a~e 



/a.j \ a — d 

3 1 2a, 



fa-[b-d+ — 



(A4) 



where we used the abbreviations 



1. Overlap molecule-substrate 



A = e Ka [J7sin(fc z a) + cos(k z a)] , 
B = c +Kb [U sin(k z b) + cos(k z b)} . 

In equation (A4) the error function erf[£] with ( e C 
arises several times. It is defined as the integral of 
the normal distribution from to £ scaled such that 
erf[±oo] = ±1: 



Let us consider the substrate case in which, for the 
sake of simplicity, we neglect in the following the spinor 
component of the substrate and atomic states. Accord- 
ing to the model given in the main text and sketched in 
Fig. 1, the substrate's wave function is given by 



erf[C] = -i f e~ t2 dt 
Jo 



and it is an entire function valid for real- and complex 
valued numbers 46 . Furthermore there holds 



(x,y,z\Sk) = -l=e +i ^ +k yy^ kz (z;z o ,0-e^ , (A5) 



& 2 

^= I e _t dt 



rf[Ca] - erf [Ci 



Ci 



Both for the wave function $j. a and the integral F kz 
the tip and the substrate cases arc obtained by the sub- 
stitutions (see also the triple- well in Fig. 1) 



sub 




tip 




where k 



x/y/z 



j?r £ x/y/z an d S is the area of the sur- 



face of the substrate on which the molecule lies. The ex- 
ponentials in Eq. (A5) stem from using no confinement 
to describe the substrate in the x and y direction and pe- 
riodic boundary conditions. Due to the large size of the 
substrate in all the three directions if compared with the 
Fermi wavelength Xp — y h 2 / (2m£p) we neglect the mo- 
mentum quantization in all three directions. By setting 
the origin of the coordinate system in R a and performing 
the Gaussian integrals in the x and y direction one easily 
obtains: 
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(Sk\a G ) =e-^I%ylSL e -^ 



Vs 



+00 



dz z * fcz (z + d, z , 0, -e )e 



=e -ifc| r fl a "Gy^ c 



x F kz (z Q ,0,-e$,aj) 



-ik,rRc 



O s (k), 



(A6) 

where fc|| • R a = k x x a + k y y a and the integral in the z 
direction has been performed with the help of Eqs. (A3) 
and (A4) . Notice the suppression of the overlap for high 
values of the parallel component of the momentum | in 
the substrate wave function given by the gaussian pref- 
actor and also the phase factor which depends on the 
position of the carbon atom R a and on kn . 

Instead of using a Gaussian p z orbital we can also use 
a Slater- type orbital 47,48 : 



(r\a s ) 



1 ( Z eS 
2^6 V a o 



(r-R a ) 



^\f-R a \ 

an 1 "i 



(A7) 

where ao = 0.53A is the Bohr radius and Z c g is a fit- 
ting parameters that takes into account the screening 
of the nuclear potential given by the core electrons. In 
Fig. 12 we show the substrate-tunnelling rates for the dif- 
ferent benzene molecular orbitals calculated according to 
Eq. (28). We compare the rates obtained using Gaussian 
and Slater-type orbitals using a distance d = 3A between 
the end of the metallic well (the substrate) and the plane 
of the molecule. As one can see the two results are in 
good agreement. The discrepancy between the two de- 
scriptions depends nevertheless on the distance d due to 
the difference in the tails of the Slater and Gaussian de- 
scriptions of the p z orbital. A good agreement is reached 
in the range of d we are interested in (d = lA — 6A). 



2. Overlap molecule-tip 

We continue with the calculation of the tip-orbital 
overlap. The atomic wave function is described again 
by the Gaussian orbitals given in Eq. (Al). The tip is 
modeled assuming a harmonic confinement in x and y 
direction, and a quantum well for the z one. The overlap 
reads: 



(x,y,z\Tk z ) 



*fc,(z;zti P ,z, 



cndi 



xe 



-(0-a;tip) 2 +(s/-s/tip) 2 ) 



(A8) 




-6 -4 -2 2 4 
AE' = AE - |Jq [eV] 

FIG. 12: (Color online) Tunneling-rates obtained by us- 
ing Slater-type orbitals (solid lines) and Gaussian orbitals 
(dashed lines). The rates are calculated for a substrate- 
molecule distance d = 3A. In the Slater-type orbital Z c g — 2. 



The overlap function is a three dimensional integral 
which, in Cartesian coordinates, reads: 

(Tk z \a G ) = n G J^J2^ 
j 

O 00 p 00 

dx Ay dze-«4(— «) 2 +(v-^) 2 +(-^) 2 ] 

X) J — OO J —OO 

x (z - d)*^ (z; z t ip, Zend, -£q) 
x e -W[( x - x iip) 2 +(y-ytip) 2 } 

(A9) 

where we have already set z a = d Va. We shift again the 
origin of the coordinates to the center of the p z orbital, 
R a , and perform the gaussian integrals in the x and y 
direction. Moreover it is convenient to introduce new 
variables describing the tip-atom distance Ax = Xti p — 
x a , Ay = ytip — Va ■ The resulting overlap function reads 



(Tk z \a G ) =n G 



E 



ft* 



a-i + 



2h 



x F kz (z tip , z Gnd , -Eq , atj) 

'■= C>T(k z , i?tip — Ra) 

(A10) 

and concludes this section dedicated to the explicit cal- 
culation of the overlap functions. 



Appendix B: The stationary density matrix 

In Eq. (38) we only gave the generic form of the station- 
ary density matrix <7 s t a t for an orbitally degenerate sub- 
space. In this section we will show how to calculate it and 
finally give the complete result for a specific example. For 
the sake of simplicity we concentrate on the transitions 
6„ O 7„, but the calculation can be easily reproduced for 
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all other transitions. The reduced density matrix <7 s t a t f° r 
the specific subspace that we are considering is composed 
of a single-element sub-block associated to the 6 particle 
ground state, a 2 x 2 sub-block associated to the sub- 
space { 1 7 g + 2 t) , 1 7 g - 2 t) } and finally a 2 x 2 sub-block 
relative to the space span{|7 9 + 2 \), \7 g — 2 J,)}. Since 
we are interested in orbital (but not spin) coherences the 
Liouvillcan is a linear operator of dimension 9x9. We 
choose the basis: 



|6 5 », 



f|7 s t; 


+2, 


+2)) 


'\7 g i 


;+2, 


f2)) 


||7 9 1; 


-2, 


-2» 


\7gl 


;-2, 


-2» 


||7 S t; 


+2, 


-2)) ' 


\7 a l 


;+2, 


-2)) 


U7 S t; 


-2, 




[\7gl 


;-2, 


f2)> 



(Bl) 



where the notation | )) denotes a vector in the density 
matrix space. We organize the tunnelling Liouvillean in 
the following form: 



stand for the diagonal elements of the tunnelling rates 
of the substrate or the tip. Moreover, the rates and the 
Fermi functions are calculated at the same energy 5E = 
E-j g — Eq 9 . The other elements of the matrix (£tun)6„7„ 
are matrices themselves. In particular: 



Cm = ^671 =fs r s (lioo 



/ T -r T ( 1 1 e + 2i ^ e~ 2i * 2 



(B3) 



are the population "rates" of the 6 particle ground state 
starting from the states \7 g l f) and \7 g l \), while 



(£tun)e 



^ ^66 ^67t ^674- 
^7t6 ^7f7t 

\ £74.6 £74,74. 



(B2) 



where = — 4 (/^r T + fg^ s ) is the depopulation rate 
of the 6 particle ground state and the coefficients r s / T 



£ 7 t6 = £716 =ftr s (110 



f£T T ( 1 1 e" 2i ^ e+ 2i " 



(B4) 



are the population "rates" of the states \7 g l f) and \7 g l i) 
starting from the state |6 ff ). Finally 



7t7t 



£ 



7171 



1 



3 + 2i</»2 



72 
/2 




1 

3 -2i0 ; 

=,+2^2 



72 
/2 



e +2i0 2 / 2 

1 




/2 e 
c 



-2i0 2 

-2102 ^2 


1 



\ 










o\ 







1 
















1 





/ 










1/ 



(B5) 



is the depopulation "rate" of the states \7 g l f) and 7 g l \) 
towards the 6 particle ground state. 

The stationary solution of the Generalized Master 
Equation Eq. (12) is found by calculating the null space 
of the Liouville operator. Here we restrict ourselves to 
the operator £t U n describing the sequential tunnelling dy- 
namics. A discussion about the relevance of the com- 
mutator with the effective Hamiltonian is left to the last 
appendix. If the leads are not superconductors, non mag- 
netic or with parallel polarization and weakly coupled to 
the molecule, the stationary density matrix is block di- 
agonal in particle number, energy and spin. Thus, the 
stationary solution which corresponds to the Liouvillean 
given in Eq. (B2) can be cast into the form: 



Cstat 




(B6) 



where the 7 particle subblocks, when written in the basis 
{|7 s +2r),|7 ff -2r)},read : 



with 



A = 
B = 



/ A Be- 2i <f> 2 

^7 9 T - <?7 9 l - I Bc +2i0 2 A 



N 

N 

r s r T (fsti-fTf£) 

N 



(B7) 



(B8) 



and the normalization N defined by the relation 
Tr(7 s tat = 1- This result is worth some further analy- 
sis. First of all it is interesting to notice that B = 
only if at least one of the following conditions is satisfied 
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i) r s = which is never happening, ii) T T = which 
holds if R tip is on the main rotational axis of benzene, 
iii) /jT / = ft / fg which is satisfied only in equilibrium 
when fix — fJ-s- This analysis shows how the interference 
between states with different angular momenta is ubiqui- 
tous in the molecular junction. Eventually, it is easy to 
prove that the eigenvalues of the stationary density ma- 
trix are <jq , A + B and A — B. The ratio between these 
eigenvalues gives a key to the physical interpretation of 
the stationary density matrix. In fact: 



where a is the distance between the carbon atoms and 
the center of the molecule, and 9 a = (27r/6)a with a = 
0, ...,5. Combining Eqs. (CI), (8), (19) and (34) we 
obtain: 



<f>l(Rti P ) = arg |^ fWxp, m P , cos(0 a - d tip )}e lWa j 

(C3) 

and, consequently: 



A + B 
A-B 

*6. 



fs 



-P(AE-ns) 



(B9) 



A-B 



which can be interpreted as follows: a^ g 
is the occupation of the 7 particle state \7 g Dr) which 
is decoupled from the tip and coupled to the 6 par- 
ticle ground states only via tunnelling events happen- 
ing at the molecule substrate interface. For this rea- 
son the ratio (Jjg/<Jd g is the same as the one obtained 
in thermal equilibrium with the substrate. On the other 
hand af g := A + B is the population of the 7 particle 
state \7 g Cr) which can exchange particles both at the 
molecule-substrate and at the molecule-tip interfaces. In 
particular, the rate of exchange for the state \7gCr) is 
double than the rate of exchange of the angular momen- 
tum states \7glr) (see Eq. (37)). The detailed balance 
gives immediately the first relation in Eq. (B9). 



Appendix C: Phase of the tunnelling amplitude 

The phase of the tunnelling amplitude between a ben- 
zene molecular orbital and a tip state plays an important 
role in the calculation of the transport characteristics of 
the STM junction. In this section we derive the approx- 
imate formula describing this phase given by Eq. (35), 
and also its limit of validity. Due to the cylindrical sym- 
metry of the tip wave function, for the overlap function 
with the atomic wave function Ox{k z , Rtip — R a ) it holds: 

OT{k Zl Rti P - Ra) = f(k z , ztip, |r t ip - r a \), (CI) 

where / is a real function (see Eq. (A10) in appendix 
A) and we have introduced cylindrical coordinates with 
the origin in the center of the molecule and the z axis 
perpendicular to the molecular plane. Every point R in 
the space is thus described by the triplet (z, r, 9) and we 
fix 6 = along the radius intersecting the atom with 
a = (see Fig. 5). Finally, we have defined r to be the 
projection of R in the plane of the molecule. It follows 
immediately that 



I rtip -r a \ = \Ja 2 + r t 2 ip - 2ar tip cos(6» tip - 9 a ), (C2) 



i(Rtip) - l9ti P = arg < ^ f(zti P , »"tip, cos 



</> Q )e 



i/0 Q 



(C4) 

where (f> a = 9 a — 9 tip ■ If now we expand / in the Taylor 
series: 



~ /(«) 

/(ztip,?"tip,COS0 Q ) = 2^ — p 



(cos cj) a ) n 



Otip,i- t ip,o) 



(C5) 

we reduce the problem to the evaluation of the functions 



9ni{9tip) = ^[cos(6» Q - 6>ti P )]"e 



iZ(8o-0ti p ) 



(C6) 



which is easily done by means of the Euler formula for 
the cosine and the binomial theorem. The solution reads: 



-i6c0 tip 



n+6a— / 
2 



[~7T 

x | cos \j^( n + 6c — I) 
x 9{n + 6c - I + 2)0(n - 6c + 1 + 2), 



(C7) 



with 9{x) = 1 if x > and elsewhere. By analyzing 
Eq. (C7) we obtain the following general properties: i) 
If tip = nir/6, with n e N, the function <? n ;(#tip) is 
real, thus Eq. (35) is exact when i? t ip is on the planes 
perpendicular to the molecule passing through the center 
of the molecule and one of the atoms or the center and 
the middle point of a carbon-carbon bond, ii) g n i = if 
n is even and g„2 = if n is odd, V0tip- iii) 3ni is real for 
n < 4 and g n i is real for n < 3. The combination of the 
observation ii) and iii) supports the validity of Eq. (35) 
on the entire space. 



Appendix D: The effective Hamiltonian 

In this section we analyze the form of the effective 
Hamiltonian H D g introduced in Eq. (39) both in the case 
of a constant interaction model or a more generic model 
for the interaction. The discussion will always be re- 
stricted to the subspace spanned by the 7 particle ground 
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states in which the effective Hamiltonian reduces to a 
4x4 matrix whose generic element is (7 ' g lr\H e fi\7 g V V) . 

Since [H m ,S z ] = it follows immediately that H e g 
is diagonal in the spin quantum number. Moreover one 
proves the following relations: 



(7 g lT\H cS \7gl'T) = (7 g lf\H cS \7 g l'f) 



(Dl) 



(7 g lT\H cS \7 g lT) = (7 g lf\H cS \7 g lf), 
which ensures i) that the two spin subblocks are iden- 



tical and ii) that the diagonal elements in each of the 
two subblocks are equal. In order to prove the relations 
given in Eq. (Dl) it is useful to introduce the symmetry 
operations C/ sp i n and U or b defined as follows: 



(D2) 



die = U orh d [a ul rh . 
The proof of the first relation in (Dl) is readily given: 



(7 g nr\H cff \7 g n'r) =-L £ [<7„nr|4r* - H m )p x (E 7g - H m )d Va \7 g n'r) 

+(7 g nr\d Va rf l ,(H In - E 7s )p x (H ni - E 7 g )d\j7 >' V) 
= ^E [{^nf\d\^UE 7g -H m )p x {E 7g - H m )d V9 \7 B n'f) 



(D3) 



U'x 

(7 g nf\d Vs Tl,{H m - E 7s ) Px (H m ~ E 7 g )d\-„\7 g n' f)] = (7 g nf\H eS \7 g n'f), 

I 



where for the second equality we have introduced the placing a — > a in the sum and remembering that Tf,, is 
7* JJ ■ 

y spin u spm 



identity operators in U spin before and after the oper- independent of the spin of the electron in the lead. The 



ators d la and d\, a . The last equality is obtained by re- second rclation in ( D1 ) is obtaincd in an analogous way: 



(7 g nr\H cH \7 g nr) =J- £ [(7„nr|4rg (E 7g - H m )p x (E 7g - H m )d lr7 \7 g 
+{7 g nT\d la Tl{H m - E 7g ) Px (H m - E 7g )dl\7 g nr) 
= ^t E [(^Ma^li^ ~ H MEj g - H m )d la \7 g 



ix 

(7 g fiT\d Trr Tl{H m - E 7g ) Px (H m - E 7g )d\j7 g , 



(D4) 



(7 g nT\H eff \7 g nT) 



where the first equality is obtained by removing the sum 
over I' since the Hamiltonian H m conserves the z projec- 
tion of the angular momentum, the second equality pro- 
ceeds instead by inserting the identities Ul lh U OI h before 
and after the operators die and d\, a . Finally, in the last 
equality, we have redefined I — >■ I and used the symmetry 
property of the rate matrices — T^. 

For the analysis of the off diagonal elements of H e s 
within a single spin subblock we have to distinguish be- 
tween the substrate and the tip case. In the substrate 
case r^, oc 6u< which directly implies that also the com- 
ponent of H e s given by the coupling to the substrate is 
diagonal and, due to the second relation in (Dl) propor- 



tional to the identity matrix and thus irrelevant for the 
dynamics of the molecule. 

Thus, let us concentrate on the tip contribution. It is 
possible to demonstrate that: 



(7 



9 -2r|^ ff |7 s 



-2r) = Ae 



-2i0 2 



Bc 



(D5) 



where we have introduced the notation H^ s to indicate 
the component of H c ff with x — T, <p\ and 4>2 are the 
phases of the tunnelling amplitudes calculated in the pre- 
vious section. Finally, A, B £ K are given by 
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A =^E [( 7 9 2 M, \rl-2(E 7g - H m )\p T (E 7g - H m )d^ 2a \7 g -2r) 

+(7 g 2r\d- 2 , \Tl_ 2 (H m - E 7g )\p T (H m - E 7g )dl\7 g -2r) , 
B =^ Rc E [( 7 9 2T \ d l \ T i3(Er g - H m )\ PT (E 7g ~ H m )d 3a \7 g -2r) 

(7 

+ (7 3 2r|d 3 . |lf 3 (iJ m - £ 73 )|p T (tf m - £Vj4a|7 fl -2r) 

I 



(D6) 



The proof of Eq. (D5) proceeds as follows. Let us start from the definition of the off diagonal matrix element: 
I 



(7 g 2 r\Hj s \7 g - 2r) =i- ]T [<7 S 2 r^r* (£ 7g - H m )p x (E 7g - H m )d Va \7 g - 2 

+ (7 9 - 2rK, CT r* (ff m - E 7g )p x (H m - E 7g )dl\7 g - 2r 
I 



(D7) 



The sums over I and /' are a priori independent and Finally, it is not difficult to prove, starting from 
run over all possible single particle angular momenta: Eq. (31), the following properties for the elements of the 
I, I' = —2, . . . , 3. The angular momentum conservation rate matrix T T : 
of H m implies, nevertheless, that the combinations which 
contribute to the sum must satisfy the condition 



2 - (-2) =1-1' (mod 6), 
which restricts the sum to the three pairs: 




(D8) 



(D9) 



rT l , = \rf l ,\e-^-M, 



(DIO) 



Combining Eq. (D7) with (D9) and (DIO), one obtains: 



{7 g 2r\H^\7 g -2r) =i- £ [i 7 9 2T \ d l K-2(E 7g - H m )\ ^^p T {E 7g - H m )d- 2lT \7 g -2 

<7 

+ (7 9 2r|cL 2(T \Tl_ 2 (H m - E 7g )\ e- 2i ^ PT (H m - E 7g )d\ a \7 g -2r) 
+ (7 9 2t\4 \Tj 3 (E 7g - H m )\ c-^p T (E 7g - H m )d 3<7 \7 g -2r) 
' (7 g 2r\d 3a \Tj 3 (H m - E 7g )\e-^ PT (H m - E 7g )d\ a \7 g -2r) 

-m- -m- \ I 1 sf~l -t , „ , , 



\ j i | i»v — ' y ' i 

-(7 3 2r|4 CT |r 3 r : _ 1 (£; 7g - J ff m )| 



l p T (E 7g -ff m )d_ lCT |7 9 -2r) 
(7 3 2r|d_ lCT \Tl^(H m - E 7g )\ e-^ PT (H m - E 7g )d\ a \7 g -2r] 



(Dll) 



from which Eq. (D5) can be easily obtained. It is now proportional to B vanishes. The eigenstates of the in- 
interesting to explore the different limits of Eq. (D5) . In teracting Hamiltonian H m coincide in fact in the con- 
the constant interaction picture, for example, the term stant interaction model with the single Slater determi- 
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nant eigenstates of the non interacting one. In practice, 
the 7 particle ground state \7 g lr) can be written as: 

416, 00), (D12) 

with 

+i 

i6 3 oo)= n n 4rio>. ( m3 ) 

l=-l r=t4 

Thus, it follows immediately that: 



(42). This choice is particularly interesting among all 
others since if #tip = nir/3 (for example when the tip is 
exactly above one of the carbon atoms) the operator L 
defined in Eq. (42) is the generator of the discrete rota- 
tions around the axis passing through the center of the 
molecule and the carbon atom closest to the tip. 

Finally let us consider under which conditions the effec- 
tive Hamiltonian commutes with the stationary density 
matrix evaluated only taking into account the tunnelling 
dynamics. By combining Eqs. (D5) and (38) one eqsily 
obtains for the 7 particle ground state with spin r: 



dl\7 g -2r)=0. 

By inserting Eq. (D14) into the second equality in (D6) 
one concludes that, in the constant interaction model, the 
effective Hamiltonian for the 7 particle ground state has 
the form: 

( K Ae^ 2 ^ 2 \ 

««>*-(J5m. k J. < D15 > 

where the hermitianicity of the i? e ff has been used. 
The constant K obtained from the direct evaluation of 
Eq. (39) is different from the off-diagonal constant A. 
Nevertheless, any contribution to the N, E, S z subblock 
of the effective Hamiltonian which is proportional to the 
unity matrix does not influence the dynamics of the sys- 
tem (see Eq. (12)). Thus we chose to set K = A which 
gives the form of the H e g given by the Eqs. (41) and 



[H eS , ffstat] = 2iB H B a sin(20 2 - ^Jtr,, (D16) 

where cr z is the third Pauli matrix and we have intro- 
duced the subscripts a and H to distinguish between the 
constants proceeding from the density matrix and the ef- 
fective Hamiltonian. In the constant interaction picture 
Bh = 0, while B a = if the tip is respecting the rota- 
tional symmetry of the molecule, i.e. i?tip is on the prin- 
cipal rotational axis of the molecule. Finally a vanish- 
ing condition can be obtained also from the phases when 
202 — 4>i = n7T - By assuming the approximate expres- 
sion for the phase given by Eq. (35) one gets 9 tl p = nir/3 
which corresponds to a tip belonging to one of the verti- 
cal mirror planes for the molecule intersecting a carbon 
atom. Notice that for these special values of 8 t i p Eq. (35) 
is exact. 

For completeness we conclude with the results regard- 
ing the 5 particle ground state. The effective Hamiltonian 
for the generic description of the Coulomb interaction 
reads: 



K 



A e -2i<Pi + Be-' 1 * 2 + Ce { ' 



K 



(H e g) 5gT - 
where A, B, C £ M. are given by 

A =^E [< 5 9 lr l4a \ T l-i( E ^ - ff m)|M#5 9 - H m )d- la \5 g -lr) 
+ (5 9 \r\d. Xa \Tl^(H m - E 5g )\p T (H m - E 7g )d\ a \7 g -lr)] , 

B =l Re Yl [(M^ML|r£>(£7 9 -H m )\ PT (E 5g -H m )d 2a \h g -\T) 

(7 

+ (5 3 lT\d 0a \Tl (H m - E 7g )\p T (H m - E 7g )4 a \5 g -lr), 
C ^ Rc E [(5 g lr\dl\T^(E 7g - H m )\ PT (E 5g ~ H m )d la \5 g -It) 



(D17) 



(D18) 



-(5 g lr\d la Irii^ - E 7g )\ PT (H m - E 7g )dl\5 g -It) 



In close analogy with the 7 particle case, one proves that B and C vanish in the constant interaction picture 



19 



and also that for 6f ip = mr/3 the effective Hamiltonian only considering the tunnelling dynamics, 
commutes with the stationary density matrix calculated 
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